Method and system for enhancing solutions to a system of linear equations

ABSTRACT

A method of significantly improving the quality of solutions to a system of linear equations. The solution to a system of linear equations is enhanced by: (1) modeling a target medium into a plurality of elements and imposing at least one localized fluctuation into the target medium; (2) measuring an output resulting from at least one localized fluctuation; and (3) processing the measured output to reconstruct a result, determining a correction filter, and applying the correction filter to the result.

This application claims the benefit under 35 U.S.C. § 120 of prior U.S.Provisional Patent Application Ser. No. 60/309,572 filed Aug. 2, 2001,entitled “A METHOD FOR FREQUENCY ENCODED SPATIAL FILTERING TO ENHANCEIMAGING QUALITY OF SCATTERING MEDIA.”

This invention was made with U.S. Government support under NIH grantnumber RO1-CA66184. The U.S. Government has certain rights in theinvention.

FIELD OF THE INVENTION

This invention relates to the field of linear equations, and moreparticularly to using correction filters to enhance solutions to asystem of linear equations such as the type of equations used in imagingof a scattering medium.

BACKGROUND

Imaging of a scattering medium relates generally to a modality forgenerating an image of the spatial distribution of properties (such asthe absorption or scattering coefficients) inside a scattering mediumthrough the introduction of energy into the medium and the detection ofthe scattered energy emerging from the medium. Systems and methods ofthis type are in contrast to projection imaging systems, such as x-ray.X-ray systems, for example, measure and image the attenuation orabsorption of energy traveling a straight line path between the x-rayenergy source and a detector, and not scattered energy. Whether energyis primarily highly scattered or primarily travels a straight line pathis a function of the wavelength of the energy and medium it is travelingthrough.

Imaging based on scattering techniques permits the use of new energywavelengths for imaging features of the human body, earth strata,atmosphere and the like that can not be imaged using projectiontechniques and wavelengths. For example, x-ray projection techniques maybe adept at imaging bone structure and other dense objects, but arerelatively ineffective at distinguishing and imaging blood oxygenationlevels. This is because the absorption coefficient of blood does notvary significantly with blood oxygenation, at x-ray wavelengths.However, infrared energy can identify the spatial variations in bloodvolume and blood oxygenation levels because the absorption coefficientat these wavelengths is a function of hemoglobin states. Otherstructures and functions can be identified by variations or changes inthe scattering coefficient of tissue exposed to infrared energy, such asmuscle tissue during contraction, and nerves during activation. Thesestructures could not be imaged by projection techniques becauseprojection techniques are not effective in measuring variations inscattering coefficients. These measures, obtainable through imagingbased on scattering techniques, such as optical tomography, haveconsiderable potential value in diagnosing a broad range of diseaseprocesses.

A typical system for imaging based on scattered energy measures,includes at least one energy source for illuminating the medium and atleast one detector for detecting emerging energy. The energy source isselected so that it is highly scattering in the medium to be imaged. Thesource directs the energy into the target scattering medium and thedetectors on the surface of the medium measure the scattered energy asit exits. Based on these measurements, a reconstructed image of theinternal properties of the medium is generated.

The reconstruction is typically carried out using “perturbationmethods.” These methods essentially compare the measurements obtainedfrom the target scattering medium to a known reference scatteringmedium; The reference medium may be a physical or a fictitious mediumwhich is selected so that it has properties that areas close as possibleto those of the medium to be imaged. Selecting a reference medium isessentially an initial guess of the properties of the target. In thefirst step of reconstruction, a “forward model” is used to predict whatthe detector readings would be for a particular source location based onthe known internal properties of the reference medium. The forward modelis based on the transport equation or its approximation, the diffusionequation, which describes the propagation of photons through ascattering medium. Next, a perturbation formulation of the transportequation is used to relate (1) the difference between the measured andpredicted detector readings from the target and reference, respectively,to (2) a difference between the unknown and known internal properties ofthe target and reference, respectively. This relationship is solved forthe unknown scattering and absorption properties of the target. Thefinal distributions of internal properties are then displayed or printedas an image.

Imaging systems and methods based on scattering techniques, such asoptical tomography systems, provide a means with which to examine andimage the internal properties of scattering media, such as theabsorption and diffusion or scattering coefficients. However, theaforementioned imaging systems and methods that recover, contrastfeatures of dense scattering media have thus far produced results havingat best modest spatial resolution. Strategies for improving imagequality are known (e.g., Newton type), but invariably these arecomputationally intensive and can be quite sensitive to initial startingconditions.

Central to the method of image formation in magnetic resonance imaging(MRI) is that there is a one-to-one correspondence between the frequencyof the measured induced current and the spatial orientation of themagnetic field gradient. Because the spatial orientation of the magneticfield gradient is known, this correspondence permits a direct assignmentof a measured response to the origin of the signal in space. In effect,the physics of the magnetic resonance phenomenon encodes a frequencysignature into the measured data that has a known spatial relationshipwith the target medium. More generally speaking, methods of this typeare known as “frequency encoded spatial filtering.”

For the foregoing reasons, there is a need for a computationallyefficient nonlinear correction method that is capable of significantlyimproving the quality of solutions to a system of linear equations suchas reconstructed images of a scattering medium.

SUMMARY OF THE INVENTION

The present invention satisfies this need by providing a method andsystem for image reconstruction and image correction that iscomputationally efficient and improves the quality of reconstructedimages of a scattering medium.

In one embodiment of the system and method of the present invention,reconstructed images of a scattering medium are enhanced by: (1)subdividing a target medium into a plurality of volume elements andassigning a modulation frequency to at least one of the volume elements'optical coefficients; (2) directing energy into the target medium fromat least one source during a period of time, and measuring energyemerging from the target medium through at least one detector; and (3)processing the measured energy emerging from the target medium toreconstruct at least one image, determining a frequency encoded spatialfilter (FESF), and applying the FESF to at least one reconstructedimage.

In another embodiment of the system and method of the present invention,a solution to a system of linear equations is enhanced by: (1) modelinga target medium into a plurality of elements and imposing at least onelocalized fluctuation into the target medium; (2) measuring an outputresulting from at least one localized fluctuation; and (3) processingthe measured output to reconstruct a result, determining a correctionfilter, and applying the correction filter to the result.

The above advantages and features are of representative embodimentsonly, and are presented only to assist in understanding the invention.It should be understood that they are not to be considered limitationson the invention as defined by the claims, or limitations on equivalentsto the claims. For instance, some of these advantages may seem mutuallycontradictory, in that they cannot be simultaneously implemented in asingle embodiment Similarly, some advantages are primarily applicable toone aspect of the invention. Thus, this summary of features andadvantages should not be considered dispositive in determiningequivalence. Additional features and advantages of the invention willbecome apparent in the following description, from the drawings, andfrom the claims.

BRIEF DESCRIPTION OF THE DRAWINGS

For a better understanding of the invention, together with the variousfeatures and advantages thereof reference should be made to thefollowing detailed description of the preferred embodiments and to theaccompanying drawings therein:

FIG. 1 is a schematic of an optical tomography system used in accordancewith the present invention;

FIG. 2 is an illustrative flowchart describing the method of the presentinvention;

FIG. 3 illustrates a target medium;

FIG. 4A illustrates the global spatial correlation between eachamplitude map and the known spatial distribution of the correspondingfrequency in the target medium, plotted as a function of the modulationfrequency (f_(m)) for the weight-transform singular-value decomposition(SVDWT) algorithm;

FIG. 4B illustrates each amplitude map's center of mass plotted as afunction of f_(m) for the SVDWT algorithm;

FIG. 5A illustrates the global spatial correlation between eachamplitude map and the known spatial distribution of the correspondingfrequency in the target medium plotted as a function of f_(m) for thecombined SVDWT with an additional matrix preconditioning operation(SVDWTWRS) algorithm;

FIG. 5B illustrates each amplitude map's center of mass plotted as afunction of f_(m) for the SVDWTWRS algorithm.

DETAILED DESCRIPTION

1. Introduction

The system and method of the present invention will be discussed inaccordance with its application to the field of optical tomography. Itis noted, however, that this methodology applies to a broad range ofproblems dealing with linear applications in which linear perturbationtheory is applied to foster a solution such as economics,quality-control, epidemiology, meteorology, or the like.

Image reconstruction methods employ computation-intensive algorithms,which are modifications of a standard linear perturbation approach toimage recovery. One of the factors that has made the development ofthese algorithms difficult in the past has been the absence of a way toquantitatively characterize the information spread function (ISF)associated with a given image reconstruction method. The term ISF usedherein refers to the precise manner in which the optical coefficientsthat actually are present in a given location of a target medium aremapped into the spatial domain of the image.

In the absence of information regarding the ISF, there is no apparentway of systematizing the process of modifying a reconstruction algorithmin response to the observed quality of its performance. In order tocharacterize the ISF for a given combination of reconstruction algorithmand reference medium the present invention utilizes the techniques foundin magnetic resonance imaging (M), which encode a frequency responseinto measurement data that has a known spatial relationship with atarget medium. Where the present invention differs from MRI is thatrather than directly applying this strategy for image formation, thepresent invention instead applies this concept to derive a frequencyencoded spatial filter (FESF) that is then applied to improve thespatial convolution of images previously recorded using other methods.

This is accomplished by recognizing that FESFs can be derived byexamination of the position-dependent temporal frequency spectraobtained from a time series of images whose optical properties in eachelement were assigned different time-varying properties. In the case ofa perfect imaging method, analysis of the time series would exactlyrecover the temporal behavior in every pixel. In practice, spatialconvolution is present, in which case the location and amplitude of theconvolving contrast feature can be determined from examination of thefrequency spectrum of the pixel data. However, by assigning temporalproperties that are uniquely distinguishable among all pixels, preciseassignment of image contrast from any one pixel to any other ispossible. The resulting information is then used as a linear operatorthat serves to rearrange (i.e., deconvolve) the contrast features of arecovered image from a test medium, thereby improving image quality.Implicit in this scheme is the assumption that the spatial convolutiondefined by the FESF is similar to the convolution present in the imageof the test medium. In principle, any number of FESFs can be derived andapplied as needed.

For illustration purposes, the present system and method is described infurther detail below with respect to an optical tomography system usedto generate images of a target scattering medium. However, it will beappreciated by those skilled in the art that the methodology of thepresent invention is applicable in image reconstruction from measureddata based on any energy source (e.g., electromagnetic, acoustic, etc.),any scattering medium (e.g., body tissues, oceans, foggy atmospheres,geological strata, and various materials, etc.), any source condition(e.g., time-independent, time-harmonic, time-resolved) and any physicalimaging domain (e.g., cross-sectional, volumetric). Accordingly, thismethodology can be extended to allow for new imaging approaches in abroad range of applications, including nondestructive testing,geophysical imaging, medical imaging, and surveillance technologies.

2. Optical Tomography System

There are many known imaging systems for collecting the measured dataused in image reconstruction in scattering media. A schematicillustration of an optical tomography system is shown in FIG. 1. Thissystem includes a computer 102, energy sources 104 and 106, a fiberswitcher 108, an imaging head 110, detectors 112, a data acquisitionboard 114, source fibers 120 and detector fibers 122.

The energy sources 104 and 106 provide optical energy, directed througha beam splitter 118, to the fiber switcher 108 and then to each of theplurality of source fibers 120 one at a time in series. The sourcefibers 120 are arranged around an imaging head 110 so that the energy isdirected into the target medium 116 at a plurality of source locationsaround the target.

The energy leaves the source fiber 120 at the imaging head 110 andenters the target medium 116 centered in the imaging head 110. Theenergy is scattered as it propagates through the target medium, emergingfrom the target medium at a plurality of locations. The emerging energyis collected by the detector fibers 122 arranged around the imaging head110. The detected energy then travels through the detector fibers 122 todetectors 112 having energy measuring devices that generate a signalcorresponding to the measurement. The data acquisition board 114receives the measurement signal for delivery to the computer 102.

This process is repeated for each source position so that a vector ofmeasures are obtained for all of the detectors and source locations. Thecomputer 102 or other suitable processing device or hardware is used toprocess the collected data and reconstruct the image as described indetail by the methods below.

3. Method

FIG. 2 is an illustrative flow chart describing the method of thepresent invention. The first step in accordance with the presentinvention is to subdivide a scattering medium for which the filterfunction will be computed into N small area or volume elements (step200). Next, a sinusoidal temporal variation is assigned to an opticalparameter (e.g., absorption and/or scattering coefficients) in eacharea/volume element, with a different frequency to each location (f₁,f₂, . . . , f_(N)) (step 210). The oscillation frequencies (i.e.,modulation frequency) in step 210 are chosen in such a way that everyratio of frequencies f_(m)/f_(n), n m, is an irrational number.

The subsequent step involves computing a time series of forward problemsolutions: I(j,t_(k)), where j=1, 2, . . . , J is the detector index andk=1, 2, . . . , K is the time-step index for the resulting dynamicmedium with N>2max(f_(m))/min(Δf_(m)) (step 220). In order to preventfrequency aliasing in step 220, the time interval between successivestates of the medium must be small relative to the reciprocal of thehighest frequency in the medium. Further, to ensure that there will besufficiently high frequency resolution in the computed time series, thetotal duration of the measurement must be long relative to thereciprocal of the smallest difference between any two assignedfrequencies. The data obtained in step 220 constitute a J×K matrix ofdetector readings.

The next step involves solving K inverse problems (i.e., reconstructingthe time series of tomographic images), one for each set of detectorreadings computed in step 220 (step 230). In step 230, each column ofthe J×K matrix of detector readings, in turn, is used to generate theleft-hand side of the equation δI=Wδx, and the corresponding x iscalculated. The data obtained in this step constitute an N×K matrix ofreconstructed optical parameters—the n n^(th) row is the time series forthe optical parameter in the n^(th) pixel or voxel.

Once the image time series is complete, the temporal discrete Fouriertransform (DFT) for each pixel of the tomographic images is computed(step 240). Subsequently, a spatial map of the DFT amplitude at eachmodulated frequency is created (step 250). The FFSF is determined byconcatenating the DFTs computed in step 240 into an array or matrix,wherein each row corresponds to the DFT amplitude in one image pixel andeach column corresponds to the DFT amplitude in one image pixel and eachcolumn corresponds to the spatial map of the amplitude at a particularfrequency (step 260). The result, in step 260, is a single linearsystem, e.g., μ_(a)*=Fμ_(a), where μ_(a) and μ_(a) are the N×1 vectorsof reconstructed and true absorption coefficients, respectively, and Fis an N×N matrix that is determined, as described subsequently, bycomparing the matrix of DFT amplitudes computed from the image timeseries to the known ideal DFT amplitude matrix. Determination ofμ_(a)*=Fμ_(a), i.e., FESF, is accomplished via a straightforward LUdecomposition (i.e., Gaussian elimination). It is noted that asingular-value decomposition (SVD) may also be used. Application of theFESF to each reconstructed image of the time series is a matter ofperforming a simple back-substitution (i.e., a spatial deconvolutioncorrection to the reconstructed images) (step 270).

3.1 Forward Model

The following discussion regarding the Forward Model (i.e., the forwardproblem) is provided to elucidate the first step of reconstruction,which is used to predict what the detector readings would be for aparticular source location based on the known internal properties of areference medium.

As discussed above, typical reconstruction techniques are based onperturbation methods that essentially relate the difference betweenpredicted detector measurements from a reference medium and detectormeasures from the target, to solve for the difference between unknownproperties of the target and known properties of the reference.Accordingly, one of the first steps in reconstruction is to select areference medium and predict the detector readings by modeling orphysical measure. Modeling the energy propagation in the scatteringmedium is done using the transport equation or its approximation, thediffusion equation. The equations describe the propagation of photonsthrough a scattering medium. For a domain having a boundary ∂Λ, this isrepresented by the expression:∇·[D(r)∇u(r)]−μ_(a)(r)u(r)=−δ(r−r _(s)),rεΛ,  (1)where u(r) is the photon intensity at position r, r_(s) is the positionof a DC point source, and D(r) and μ_(a)(r) are the position-dependentdiffusion and absorption coefficients, respectively. Here the diffusioncoefficient is defined as D(r)=1/{3[μ_(a)(r)+μ_(s) (r)}, where μ_(s) (r)is the reduced scattering coefficient. Using this equation, the energyemerging from the reference medium at each detector location for eachsource location is predicted. The transport or diffusion equations arealso the basis for formulating the perturbation or inverse formulationused in reconstruction.3.2 The Inverse Formulation

The following discussion regarding the Inverse Formulation (i.e., theinverse problem), is provided to elucidate the second step ofreconstructing the time series of tomographic images.

As discussed above, reconstruction of a cross-sectional image of theabsorption and/or scattering properties of the target medium is based onthe solution of a perturbation or inverse formulation of the radiationtransport or diffusion equation. The perturbation method assumes thatthe composition of the unknown target medium deviates only by a smallamount from a known reference medium. This reduces a highly non-linearproblem to one that is linear with respect to the difference inabsorption and scattering properties between the target medium underinvestigation and the reference medium. The resulting optical inverse orperturbation formulation is based on the normalized difference methodand has the following form:W _(r) ^((μ) ^()·δμ) _(a) +W _(r) ^((D)) ·δD=δμ _(r),  (2)where δμ_(a) and δD are the vectors of cross-sectional differencesbetween the optical properties (absorption and diffusion coefficients,respectively) of a target (measured) medium and of a reference (computedor measured) medium used to generate the initial guess; W_(r) ^((μ) ^(a)⁾ and W_(r) ^((D)) are the weight matrices describing the influence thatlocalized perturbations in the absorption and diffusion coefficients,respectively, of the selected reference medium have on the surfacedetectors; and μu_(r) represents a normalized difference between twosets of detector readings, which is defined by the equation:$\begin{matrix}{{\left( {\delta\quad u_{r}} \right)_{i} = {{\left( \frac{\left( u_{1} \right)_{i} - \left( u_{2} \right)_{i}}{\left( u_{2} \right)_{i}} \right)\left( u_{r} \right)_{i^{*\quad}}\quad i} = 1}},2,\ldots\quad,{M.}} & (3)\end{matrix}$Here, u_(r) is the computed detector readings corresponding to theselected reference medium, u₂ and u₁ represent two sets of measured data(e.g., background vs. target, time-averaged mean vs. a specific timepoint, etc.) and M is the number of source-detector pairs in the set ofmeasurements.3.3 Weight Matrix Scaling

The following discussion regarding Weight Matrix Scaling is provided toelucidate the scaling of the weight matrices arrived at in the inverseproblem.

The effect of scaling the weight matrix is to make it more uniform,which can often serve to improve its conditioning. A scaling approachthat scales each column of W_(r) ^((μ) ^(a) ⁾ and W_(r) ^((D)) to theaverage value of the column vector is used. However, it should beunderstood that any of the known scaling approaches could be adopted.The form of the resulting new weight matrices is:{tilde over (W)} _(r) ^((k)) =W _(r) ^((k)) ·R ^((k)),  (4)where k can be μ_(a) or D, and R^((k)) is the normalizing matrix whoseentries are: $\begin{matrix}{\left( R^{(k)} \right)_{ij} = \left\{ {{\begin{matrix}\frac{1}{\frac{1}{M}{\sum\limits_{m = 1}^{M}\left( W_{r}^{(k)} \right)_{mj}}} & {{j = i},} \\0 & {{j \neq i},}\end{matrix}\quad i},{j = 1},2,\ldots\quad,N,} \right.} & (5)\end{matrix}$in which N is the number of elements used in discretizing the domain Λ.The resulting system equation is:{tilde over (W)}_(r) ^((μ) ^(a) ⁾·δ{tilde over (μ)}_(a) {tilde over (W)}_(r) ^((D))·δ{tilde over (D)}=δμ_(r),  (6)

-   -   where δ{tilde over (μ)}_(a)=[R^((μ) ^(a) ^()]) ⁻¹·δμ_(a) and        δ{tilde over (D)}=[R^((D))]⁻¹·δD. Note that R^((k)) is a        diagonal matrix (Eq. 5) all of whom main diagonal elements are        non-zero (Eq. 5); consequently it has a well-defined inverse,        the computation of which is a trivial matter.        4. Frequency Encoded Spatial Filter (FESF)

The following discussion regarding FESF is provided to elucidatedetermination of the FESF, which is used to reconstruct correctedimages.

The amplitude “spatial maps” produced by image reconstruction and thecomputations of the DFTs in actuality are strings of numbers, each beingthe amplitude, at one particular frequency, assigned by thereconstruction algorithm to one of the FEM nodes (i.e., the number orvertices, or points where three or more elements come together). Theentire set of amplitude maps can be concatenated into a matrix:${A_{i} = \begin{bmatrix}A_{11} & A_{12} & \ldots & A_{1N_{n}} \\A_{21} & A_{22} & \ldots & A_{2N_{n}} \\\vdots & \vdots & ⋰ & \vdots \\A_{N_{f}1} & A_{N_{f}2} & \ldots & A_{N_{f}N_{n}}\end{bmatrix}},$where N_(f) is the number of frequencies (=number of finite elements)and N_(n) is the number of FEM nodes. This is important in understandingthe determination of the FESF, because in practice the number of nodesinvariably is smaller than the number of elements. Then A_(i) (i for“mages”) is not a square matrix, but has approximately half as many rowsas columns.

As such, a second matrix can be written, which tells us exactly whereeach modulation frequency actually was present in the medium:$B_{m} = {\begin{bmatrix}B_{11} & B_{12} & \ldots & B_{1N_{n}} \\B_{21} & B_{22} & \ldots & B_{2N_{n}} \\\vdots & \vdots & ⋰ & \vdots \\B_{N_{f}1} & B_{N_{f}2} & \ldots & B_{N_{f}N_{n}}\end{bmatrix}.}$The matrix B_(m) must be sparse (i.e., most of its elements are zeroes),because each f_(m) is assigned to only one of the medium's finiteelements. In fact, every row of B_(m), which contains hundreds orthousands of elements altogether (i.e., N_(n)=10³−10⁴), has exactlythree (in the case of two-dimensional media) or four (three-dimensionalmedia) elements that are not zero. The number is three or four due tothe use of triangular or tetrahedral elements, so each element isbounded by three or four nodes.

If the previously described reconstruction process were perfect,A_(i)=B_(m). This ideal result, however, is not achieved in practice.Thus one must make some type of assumption regarding the nature of thefunction that transforms A_(i) into B_(m). The one made in the presentinvention is that the frequency spectrum present at any one node (i.e.,pixel) in the images is a linear function of the frequencies present atall nodes in the medium. Mathematically, this is stated by:TA_(i)=B_(m), where T is a N_(n)×N_(n) (i.e., square) matrix. Inpractice, one wants the transformation to go in the other direction,that is, starting from B_(m), produce something that is as close aspossible to the true A_(i). Thus, computation of the filter that willactually be used in practice is accomplished by solving the matrixequation A_(i)=FB_(m), where F also is a N_(n)×N_(n) matrix. T and F areinverses of each other.

It is noted that the total number of elements in F is smaller than thenumber in B_(m), because N_(n)<N_(f). This means that perfect correctionwill not occur when applying this method, because there aren't enoughcorrection terms to go around. This occurs, not because of the assumedlinear relation between A_(i) and B_(m), but because there is anunavoidable loss of information associated with mapping the frequenciesin N_(f) elements into the smaller number, N_(n), of nodes.

The FESF that is computed in this way has quality-control utility as away of quantifying the accuracy of reconstruction algorithms. The FESFmay further be used as an image enhancing tool if it is employed inconjunction with data obtained from different experimental media fromthe one used to generate the filter. In this scenario, uponreconstruction of a set of images I₁, I₂, etc., then, to the extent thatthe filter function is not strongly dependent on the medium'sproperties, the spatial accuracy of the reconstruction can be improvedby computing FI₁, FI₂, etc.

5. Demonstration Results

The following example is presented to illustrate features andcharacteristics of the present invention, and is provided solely toassist in explanation of a demonstration of the invention and is notintended to be construed as limited thereto.

A demonstration of the utility of the FESF is described in the foregoingexample. As discussed below the FESF has been applied to two differentimage time series, both obtained from the same sets of detector readingsbut employing different varieties of a reconstruction algorithm. Inprinciple the reconstruction methods employed should produce identicalresults since there is no self-evident a priori reason for choosing touse one rather than the other. However, application of the FESF methodindicates that one variety of reconstruction methods can producespatially accurate images of perturbations at any location of themodeled medium, and the other can not do so. Accordingly, the computedISF for either algorithm affords a way of applying a spatialdeconvolution correction to a reconstructed image.

FIG. 3 illustrates a regularly-shaped two-dimensional medium (i.e., thetarget medium). As shown in FIG. 3, the medium is a homogeneous disk of8 cm diameter, with optical coefficient values of _(a)=0.06 cm⁻¹,_(s)=Dcm⁻¹. For more convenient solution of the forward and inverseproblems, the mathematical boundary of the disk was extended 0.5 cmbeyond that of the “physical” medium, as indicated in FIG. 3. Thecoefficient values in the extended region were the same as those of the“physical” medium. Sixteen equally spaced, unit-strength, homogeneouspoint sources were placed in the medium at the indicated positions onthe physical boundary.

The numbers of finite elements and nodes in the indicated mesh are 1604and 850, respectively, and the smallest and largest element areas are0.026 cm² and 0.073 cm² (mean±standard deviation=0.040±0.006 cm²).Sinusoidal modulation was imposed on the absorption coefficient in eachelement. A unique f_(m) was assigned to each, while the amplitude waseverywhere 0.006 cm⁻¹ (i.e., 10% of the mean value). For thispreliminary study, the elements' scattering coefficients were notmodulated in time.

To ensure that the resolution bandwidth was smaller than the smallestdifference between f_(m)s and the Nyquist frequency was greater than thelargest f_(m), a time series of ten thousand sets of tomographicdetector readings was computed, with t=0.01 s. Image reconstruction wascarried out with two algorithms, both based on an SVD of the imageoperator matrix. The first algorithm used was the previously describedweight-transform SVDWT method. The second reconstructionmethod—SVDWTWRS—combined SVDWT with an additional matrix preconditioningoperation, in which each equation was scaled so that all rows of theweight matrix had the same sum.

Two types of analysis were performed on the 1,604 DFT amplitude mapsproduced during image reconstruction. First, the global spatialcorrelation was computed between each amplitude map and the knownspatial distribution of the corresponding frequency in the target medium(ideal result: correlation exactly equal to 1.0 at all frequencies).Second, the coordinates of each map's center-of-mass were computed, fromwhich we easily determined its displacement from the geometric centroidof the finite element whose _(a) was modulated at the correspondingfrequency (ideal result: displacement exactly equal to 0.0 at allfrequencies). These two quantities are plotted, as a function of f_(m)(or, equivalently, location in the target medium), for the SVDWTalgorithm in FIGS. 4A and 4B, and for the SVDWTWRS algorithm in FIGS. 5Aand 5B. The lighter-colored curve in FIGS. 4A and 4B, and 5A and 5B arederived from the unfiltered FT amplitude spatial distributions. Thedarker curves are the results obtained when the calculations wererepeated after we made the best-possible correction consistent with thetheoretical model described above, according to which the amplitude mapsderived from the reconstructed images are a simple linear transformationof the true spatial distributions present in the target medium.

Inspection of FIGS. 4A and 4B, and 5A and 5B reveals that each plottedfunction exhibits a qualitative change in behavior after the 400^(th)f_(m). The change is simply a consequence of the fact that the first 400finite elements all were located in the zone (see FIG. 3) lying betweenthe physical and extended boundaries, i.e., outside the ring of sourcesand detectors. Closer inspection of FIGS. 4A and 4B reveals that bothspatial accuracy measures fall particularly far from their ideal valuesfor those finite elements corresponding to roughly the 800^(th) through1100^(th) f_(m). These elements are the ones that lay in the centralregion of the target medium. That is, the SVDWT algorithm reconstructedimages that were strongly distorted spatially, with the absorptioncoefficient values of the central region significantly displaced towardthe surface while those of the more peripheral region were recoveredwith considerably greater accuracy. In contrast, the spatial correlationand centroid displacement are considerably more spatially uniform forthe amplitude maps derived from the images reconstructed by the(preconditioned) SVDWTWRS algorithm. This is a significant observation,as the two reconstruction variants theoretically should yield the samesolution when both operate on a given set of detector data. Finally, itis seen that in each panel of FIGS. 4A and 4B, 5A and 5B, most points onthe dark (corrected images) curve lie closer to the ideal value thanthose on the light (uncorrected mages) curve. This demonstrates thepossibility that information in the ISF could be used to performpost-reconstruction enhancement of the images' spatial accuracy.

It should be understood that the above description is onlyrepresentative of illustrative embodiments. For the convenience of thereader, the above description has focused on a representative sample ofpossible embodiments, a sample that is illustrative of the principles ofthe present invention. The description has not attempted to exhaustivelyenumerate all possible variations. That alternate embodiments may nothave been presented for a specific portion of the invention, or thatfurther undescribed alternate embodiments may be available for aportion, is not to be considered a disclaimer of those alternateembodiments. Other applications and embodiments can be conceived bythose without departing from the spirit and scope of the presentinvention. It is therefore intended, that the invention is not to belimited to the disclosed embodiments but is to be defined in accordancewith the claims that follow. It can be appreciated that many of thoseundescribed embodiments are within the scope of the following claims,and others are equivalent.

1. A method of enhancing reconstructed images of a scattering medium,comprising: subdividing a first target medium into a plurality of volumeelements; assigning a modulation frequency to at least one of the volumeelements' optical coefficients; directing energy into the first targetmedium from at least one source during a period of time; measuringenergy emerging from the first target medium through at least onedetector; processing the measured energy emerging from the first targetmedium to reconstruct at least one image; determining a frequencyencoded spatial filter (FESF); and applying the FESF to at least onereconstructed image of the first target medium.
 2. The method accordingto claim 1, wherein the modulation frequency is assigned to anabsorption coefficient of a volume element.
 3. The method according toclaim 1, wherein the modulation frequency is assigned to a scatteringcoefficient of a volume element.
 4. The method according to claim 1,wherein the measured energy emerging from the first target medium isprocessed by employing a perturbation method.
 5. The method according toclaim 4, wherein the perturbation method employed uses a forward problemsolution to reconstruct the tomographic images.
 6. The method accordingto claim 4, wherein the perturbation method employed uses the inverseproblem to reconstruct the tomographic images.
 7. The method accordingto claim 1, wherein determining the FESF, comprises: computing thetemporal discrete Fourier transform of the reconstructed tomographicimages; and processing the computed temporal discrete Fourier transformto determine the amplitude at a modulation frequency associated with itscorresponding volume element.
 8. The method according to claim 1,wherein the FESF is applied to at least one reconstructed image byperforming a simple matrix multiplication.
 9. The method according toclaim 1, further comprising: directing energy into a second targetmedium from at least one source during a period of time; measuringenergy emerging from the second target medium through at least onedetector; processing the measured energy emerging from the second targetmedium to reconstruct at least one image; and applying the FESFdetermined from the first target medium to at least one reconstructedimage of the second target medium.
 10. A system for enhancingreconstructed images of a scattering medium, comprising: means forsubdividing a first target medium into a plurality of volume elements;means for assigning a modulation frequency to at least one of the volumeelements' optical coefficients; means for directing energy into thefirst target medium from at least one source during a period of time;means for measuring energy emerging from the first target medium throughat least one detector; means for processing the measured energy emergingfrom the first target medium to reconstruct at least one image; meansfor determining an FESF; and means for applying the FESF to at least onereconstructed image of the first target medium.
 11. The method accordingto claim 10, wherein the modulation frequency is assigned to anabsorption coefficient of a volume element.
 12. The method according toclaim 10, wherein the modulation frequency is assigned to a scatteringcoefficient of a volume element.
 13. The method according to claim 10,wherein the measured energy emerging from the first target medium isprocessed by employing a perturbation method.
 14. The method accordingto claim 13, wherein the perturbation method employed uses a forwardproblem solution to reconstruct the tomographic images.
 15. The methodaccording to claim 13, wherein the perturbation method employed uses theinverse problem to reconstruct the tomographic images.
 16. The methodaccording to claim 10, wherein determining the FESF, comprises: meansfor computing the temporal discrete Fourier transform of thereconstructed tomographic images; and means for processing the computedtemporal discrete Fourier transform to determine the amplitude at amodulation frequency associated with its corresponding volume element.17. The method according to claim 10, wherein the FESF is applied to atleast one reconstructed image by performing a simple matrixmultiplication.
 18. The method according to claim 10, furthercomprising: means for directing energy into a second target medium fromat least one source during a period of time; means for measuring energyemerging from the second target medium through at least one detector;means for processing the measured energy emerging from the second targetmedium to reconstruct at least one image; and means for applying theFESF determined from the first target medium to at least onereconstructed image of the second target medium.
 19. A program stored ona computer readable medium and executable by a processor, comprising:instruction code which, when executed by the processor subdivides afirst target medium into a plurality of volume elements; instructioncode which, when executed by the processor assigns a modulationfrequency to at least one of the volume elements' optical coefficients;instruction code which, when executed by the processor directs energyinto the first target medium from at least one source during a period oftime; instruction code which, when executed by the processor measuresenergy emerging from the first target medium through at least onedetector; instruction code which, when executed by the processorprocesses the measured energy emerging from the first target medium toreconstruct at least one image; instruction code which, when executed bythe processor determines an FESF; and instruction code which, whenexecuted by the processor applies the FESF to at least one reconstructedimage of the first target medium.
 20. A program, according to claim 19,further comprising: instruction code which, when executed by theprocessor directs energy into a second target medium from at least onesource during a period of time; instruction code which, when executed bythe processor measures energy emerging from the second target mediumthrough at least one detector; instruction code which, when executed bythe processor processes the measured energy emerging from the secondtarget medium to reconstruct at least one image; and instruction codewhich, when executed by the processor applies the FESF determined fromthe first target medium to at least one reconstructed image of thesecond target medium.
 21. A method of enhancing the solution to a systemof linear equations, comprising: modeling a first target medium into aplurality of elements; imposing at least one localized fluctuation tothe target medium; measuring an output resulting from at least onelocalized fluctuation; processing the measured output to reconstruct aresult; determining a correction filter; and applying the correctionfilter to the result.